Wang et al. Nanoscale Research Letters 2014, 9:300 
http://www.nanoscalereslett.eom/content/9/1/300 



O Nanoscale Research Letters 

a SpringerOpen Journal 



NANO EXPRESS Open Access 



Transient viscoelasticity study of tobacco mosaic 
virus/Ba 2+ superlattice 

Haoran Wang 1 , Xinnan Wang 1 ", Tao Li 2 and Byeongdu Lee 2 



Abstract 

Recently, we reported a new method to synthesize the rod-like tobacco mosaic virus (TMV) superlattice. To explore 
its potentials in nanolattice templating and tissue scaffolding, this work focused the viscoelasticity of the superlattice 
with a novel transient method via atomic force microscopy (AFM). For measuring viscoelasticity, in contrast to 
previous methods that assessed the oscillating response, the method proposed in this work enabled us to 
determine the transient response (creep or relaxation) of micro/nanobiomaterials. The mathematical model and 
numerical process were elaborated to extract the viscoelastic properties from the indentation data. The adhesion 
between the AFM tip and the sample was included in the indentation model. Through the functional equation 
method, the elastic solution for the indentation model was extended to the viscoelastic solution so that the time 
dependent force vs. displacement relation could be attained. To simplify the solving of the differential equation, 
a standard solid model was modified to obtain the elastic and viscoelastic components of the sample. The 
viscoelastic responses with different mechanical stimuli and the dynamic properties were also investigated. 
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Background 

The recognition of tobacco mosaic virus (TMV) since the 
end of nineteenth century [1] has sparked innumerable re- 
search towards its potential applications in biomedicine 
[2,3] and biotemplates for novel nanomaterial syntheses 
[4,5]. A TMV is composed of a single-strand RNA that is 
coated with 2,130 protein molecules, forming a special 
tubular structure with a length of 300 nm, an inner diam- 
eter of 4 nm, and an outer diameter of 18 nm [6]. The 
TMVs observed under a microscope can reach several 
tens of microns in length due to its unique feature of 
head- to-tail self-assembly [7]. Practically useful properties 
of the TMVs include the ease of culture and broad range 
of thermal stability [8]. Biochemical studies have shown 
that the TMV mutant can function as extracellular matrix 
proteins, which guide the cell adhesion and spreading [8]. 
It has also been confirmed that stem cell differentiation 
can be enhanced by both native and chemically modified 
TMV through regulating the gene's expression [9-11]. 
Moreover, TMV can be electrospun with polyvinyl alcohol 
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(PVA) into continuous TMV/PVA composite nanofiber to 
form a biodegradable nonwoven fibrous mat as an extra- 
cellular matrix mimetic [12]. 

Very recently, we have reported that the newly synthe- 
sized hexagonally packed TMV/Ba 2+ superlattice mater- 
ial can be formed in aqueous solution [13,14]. Figure 1 
shows the schematic of the superlattice formation by 
hexagonal packing of TMVs, triggered by Ba ions, and 
the images observed from field emission scanning elec- 
tron microscopy (FESEM) and atomic force microscopy 
(AFM). The sample we used for this experiment was 
tens of microns in length, 2-3 microns in width (from 
FESEM), and several hundred nanometers in height 
(from AFM height image). It is known that the superlat- 
tice exhibits physical and mechanical properties that dif- 
fer significantly from its constituent materials [15-20]. 
The study on the viscoelastic properties of the TMV- 
derived nanostructured materials is still lacking despite 
the availability of the elastic property of the TMV and 
TMV-based nanotube composites [7]. The viscoelasticity 
of micro/nanobioarchitecture significantly affects the tis- 
sue regeneration [21] and repair [22], cell growth and 
aging [23], and human stem cell differentiation [24] as 
well as the appropriate biological functions of the 
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membranes within a specific nanoenvironment [25]; in 
particular, the viscoelasticity of some viruses plays key 
roles in the capsid expansion for releasing nucleic acid 
and modifying protein cages for vaccine delivery pur- 
poses [26]. Specifically, for TMV superlattice, its nano- 
tube structure makes it a perfect biotemplate for 
synthesizing nanolattices that have been confirmed to 
possess extraordinary mechanical features with ultralow 
density [27,28]. Considering the biochemical functions 
of the TMV, its superlattice is an excellent candidate for 
bone scaffolding where the time-dependent mechanical 
properties become determinant [29], and research on 
scaffolding materials remains a hotspot [30]. Apart from 
contributing to the application of TMV superlattice, this 
work also pioneered in the viscoelasticity study of virus 
and virus-based materials. By far, most literature on viral 
viscoelasticity has been focused on the dynamic proper- 
ties of virus suspensions or solutions [31-34]. One of the 
rare viscoelasticity studies on individual virus particle is 
the qualitative characterization of the viscoelasticity of the 
cowpea chlorotic mottle virus [26] using quartz crystal 
microbalance with dissipation technique, which presents 
only the relative rigidity between two samples. To date, lit- 
tle literature is available on the quantitative study of the 
viscoelasticity of individual virus/virus-based particles. 
Considering the potential uses of TMV/Ba 2+ superlattice, 
its viscoelastic properties and responses under different 
mechanical stimuli need to be investigated. 

A number of techniques for measuring the viscoelasticity 
of macro-scale materials have been used. A comprehensive 



review of those methods can be found in the literature [35] 
that addresses the principles of viscoelasticity and experi- 
mental setup for time- and frequency-domain measure- 
ments. When the sample under investigation is in micro or 
even nanometer scale, however, the viscoelastic measure- 
ments become much more complicated. In dynamic 
methods, shear modulation spectroscopy [36] and mag- 
netic bead manipulation [37] are two common methodolo- 
gies to obtain the micro/nanoviscoelastic properties. To 
improve the measurement accuracy, efforts have been 
made to assess the viscoelasticity of micro/nanomaterials 
using contact-resonance AFM [38-41]. The adhesion be- 
tween the AFM probe tip and sample, however, is usually 
neglected. Furthermore, in order for the dynamic method 
to obtain a sinusoidal stress response, the applied strain 
amplitude must be kept reasonably small to avoid chaotic 
stress response and transient changes in material proper- 
ties [42]. In addition, the dynamic properties are frequency 
dependent, which is time consuming to map the viscoelas- 
ticity over a wide range of frequencies. An alternative way 
to measure the viscoelastic response of a material is the 
transient method. Transient indentation with an indenter 
was developed based on the functional equation methods 
[43], where the loading or traveling histories of the in- 
denter need to be precisely programmed. 

In this study, the viscoelastic properties of the TMV/ 
Ba 2+ superlattice were investigated using AFM-based 
nanoindentation. AFM has the precision in both force 
sensing and displacement sensing, although it lacks the 
programing capability in continuous control of force and 
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displacement. To realize the transient indentation in AFM, 
we introduced a novel experimental method. Viscoelas- 
tic nanoindentation theories were then developed based 
on the functional equation method [44]. The adhesion 
between the AFM tip and the sample, which signifi- 
cantly affected the determination of the viscoelastic 
properties [45], was included in the indentation model 
[20]. The viscoelastic responses of the sample with re- 
spect to different mechanical stimuli, including stress 
relaxation and strain creep, were further studied. The 
transition from transient properties to dynamic proper- 
ties was also addressed. 

Methods 

The TMV/Ba 2+ superlattice solution was obtained from 
the mixture of the TMV and BaCl 2 solution (molar ratio 
of Ba 2+ /TMV = 9.2 x 10 4 :1) as stated in the reference 
[13]. It was further diluted with deionized water (volume 
ratio 1:1). A 10-|iL drop of the diluted solution on a sili- 
con wafer was spun at 800 rpm for 10 s to form a 
mono-layer dispersion of the sample. The sample was 
dried for 30 min under ambient conditions (40% R.H., 
21°C) for AFM (Dimension 3100, Bruker, Santa Barbara, 
CA, USA) observation and subsequent indentation tests. 

The sample was observed with FESEM and AFM. The 
indentation was performed using the AFM nanoindenta- 
tion mode (AFM probe type: Tapl50-G, NanoAndMore 
USA, Lady's Island, SC, USA). The geometry of the can- 
tilever was precisely measured using FESEM (S-4700, 
Hitachi, Troy, MI, USA), with a length of 125 [im, width 
of 25 (im, and thickness of 2.1 (im. To accurately meas- 
ure the tip radius, the tip was scanned on the standard 
AFM tip characterizer (SOCS/W2, Bruker) and the 
scanned data was curve fitted using PSI-Plot (Poly Software 
International, Orangetown, NY, USA). The tip radius 
calculated to be 12 nm. For a typical indentation test, 
the tip was pressed onto the top surface of the sample 
until a predefined force of -100 nN. The cantilever end 
remained unchanged in position during the controlled 
delay time. A series of indentations of the same prede- 
fined indentation force and different delay times were 
performed to track the viscoelastic responses. A 10-min 
time interval of the two consecutive indentations was 
set for the sample to fully recover prior to the next in- 
dentation. The sample drift was minimized by turning 
off the light bulb in the AFM controller during scanning 
to keep the AFM chamber temperature constant and by 
shrinking the scan area gradually down to 1 nm x 1 nm 
on the top surface of the sample to rid the scanner piezo 
of the hysteresis effect. 

Mathematical formulation 

Derived from the functional equation method and the 
standard solid model (shown in the Appendix'), the 



differential equation governing the contact behavior of 
viscoelastic bodies can be obtained as 

^Ai^\[F(t) + 2nwR] (1) 




where F(t) is the contact force history, S(t) is the inden- 
tation depth history, R is the nominal radius of the two 
contact spheres, w is the adhesive energy density, A t and 
Bi (i = 0, 1, 2) are the parameters determined by the 
mechanical properties of two contact bodies, and the 
calculation of all these parameters can be found in the 
Appendix/ 

The elastic moduli E 1 and E 2 and viscosity rj in Figure 2 
are implicitly included in the above differential equation. 
To determine E\, E 2 , and tj, besides experimental data 
for t and F, the function of the force history F(t) is also 
required. The experimental data of t and F can be ob- 
tained as indicated in Figure 3. The force relaxation can 
be found in Figure 3a where the force decrease between 
the right ends of extension and retraction curves. By 
mapping the force decrease at different delay times as 
shown using the red asterisks in Figure 3b, the force re- 
laxation curve can be obtained, which decreases from 
104 to 40 nN. The function of F(t) can be obtained from 
Equation (1). Not only is Equation (1) applicable for the 
standard solid model in Figure 2(a) where it is derived 
from, but also it can be used for the modified standard 
solid model in Figure 2(b) where the elastic component 
of Ei is replaced by two elastic components in series. 
With this modification, the deflection of the cantilever 
can be incorporated into the deformation of the imagin- 
ary sample which is represented by the modified stand- 
ard solid model where the elastic component of E lc in 
Figure 2(b) denotes the cantilever and the rest compo- 
nents denote the TMV/Ba 2+ superlattice. 

During each indentation, the vertical distance between 
the substrate and the end of the cantilever remains con- 
stant. Therefore, as the sample deformation or the in- 
dentation depth increases, the corresponding cantilever 
deflection Ad or the normal indentation force decreases. 
During this process, the force on the system decreases 
while the sample deformation S increases to compensate 
the decreased cantilever deflection. Therefore, the 
change of the cantilever deflection is equal to change of 
the sample deformation during indentation, as is shown 
in Figure 4. As such, S in Equation (1) represents the 
relative approach between the cantilever end and the 
substrate, which incorporates the deformation of both 
the sample and the cantilever. 
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To be clearer, S is substituted by D which represents 
the combined deformation. The relative approach, D, 
can be written as 



D(t) = D 0 H(t) 



(2) 



where H(t) is the Heaviside unit step function and D 0 is 
the relative approach between the substrate and the end 
of the cantilever. 
Thus, Equation (1) can be rewritten as 



5> 



{F{t) + 2nwR) = (WR/3)4^E B '3?j /f W 

(3) 



Applying Laplace transform, it yields 

(Ao+A 1 s+A 2 ^)(f(s)+^) 
4V^/3)d|(£ 0 + Bis + B 2 s 2 ) - 



(4) 



where a function with <A ' denotes Laplace-transformed 
function in s domain. 

Performing inverse Laplace transform, the viscoelastic 
equation of AFM-based indentation becomes 



F(t) 



4 \ v o R 



(A r e- at +B r e-' it + C r )-2nwR 



(5) 
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Figure 3 Indentation force, (a) Indentation force decrease with delay time set as 100 ms, 200 ms, 500 ms, and 1,000 ms, respectively, (b) 
Indentation force vs. time data from experiment measurement and fitted curve from the indentation equation. 
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where 



T1G\K\ 



(3K X + 4G 1 )(3/C 1 G 1 + 3KiG 2 + 4GiG 2 ) 
3GiG 2 



Gi + G 2 V Gi + 3/Ci G 2 + 4Gi G 2 y ' 

= Gi + G 2 G 2 3/CiGi 

rj 1 1 rj //(3/C 1 +4G 1 ) 

Solution to AFM-based indentation equation 

It is observed from Figure 3 that the initial indentation 
force at t = 0 was measured to be 104.21 nN, then the 
force started to decrease and then remained constant at 
38 nN after ~ 5,000 ms. The force decrease shown as red 
asterisks in Figure 3b fits qualatitatively well with the ex- 
ponential function of Equation (5). E 1} E 2 , and rj, corre- 
sponding to the mechanical property parameters in 
Figure 2(a), can then be determined by fitting Equation 
(5) with the experimental data. 

From the indentation data, D 0 is obtained to be 
78.457 nm. The pull-off force, 2nwR, calculated by aver- 
aging the pull-off forces of multiple indentations on the 
sample, is 16 nN. In comparison with the radius of the 
AFM tip, the surface of the sample can be treated as a 
flat plane. Hence, the nominal radius R = R tip = 12 nm. 

By invoking the force values at t = 0, £=<*>, and any 
intermediate point into Equation (5), the elasticity and 
viscosity components can be determined to be E x = 
32.0 MPa, E 2 = 21.3 MPa, and rj = 12.4 GPa ms. The co- 
efficient of determination R 2 of the viscoelastic equation 
and the experimental data is -0.9639. 

Since the stress relaxation process is achieved by mod- 
eling a combination of the cantilever and the sample, 



the viscoelasticity of the sample can be obtained by sub- 
tracting the component of the cantilever from the re- 
sults. The cantilever, acting as a spring, is in series with 
the sample, represented by a standard solid model. The 
schematic of the series organization is shown in Figure 2 
(b). Thus the component of E x comprises of E Xs repre- 
senting the elastic part from the sample and E lc repre- 
senting the elastic part from the cantilever. To clarify 
the sources of the components in the modified standard 
solid model, E 2 , v 2 , and rj in Figure 2(a) are now respect- 
ively denoted by E 2s , v 2s) and rj s in Figure 2(b), where the 
subscript V denotes the sample. 

At the onset of indentation, only the spring with elastic 
modulus of Ei takes the instantaneous step load; therefore, 
the elastic modulus of E ls can be determined from the ex- 
perimental data of zero-duration indentation. Applying 
the DMT model [46] with the force- displacement rela- 
tionship of the cantilever, 



E — ^ ^cantilever 



(6) 



we can obtain the elastic equation of AFM-based 
indentation 



_ F (F + 2ttwR\* 



(7) 



where k is the spring constant of the cantilever, which is 
5 nN/nm based on Sader's method [47] to calibrate /c, 
^cantilever is the cantilever deflection, and 8 is recorded 
directly as the Z-piezo displacement by AFM. 

The elastic modulus of E ls can be calculated by fitting 
the DMT-model-based indentation equation with ex- 
perimental data as shown in Figure 5. For simplicity, 
modification was done to the indentation equation and 
the experimental data, whose details can be found 
in reference [20]. The fitted elastic modulus of E ls 
is -2.14 GPa with a coefficient of determination of 0.9948. 
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Figure 5 Indentation force data as a function of Z-piezo 
displacement, a comparison of experimental measurement and 
fitted results. 



Results and discussion 

Based on the solution obtained, the viscoelastic equation 
of AFM-based indentation for TMV/Ba 2+ superlattice is 
written as 



F(t) = 3.2098 ( 0.0076"^ + 0.0136e"^ + 0.0168) 



-16 

(8) 



The force decrease curve is shown in Figure 3b with 
the experimental data. 

Specifically, for the TMV/Ba 2+ superlattice whose 
viscoelastic behavior is simulated by a standard solid 
model, the differential equation governs its stress-strain 
behavior and becomes 



E 2s Ei s E 2s / \ • 

o H a = £+ [Eis + E 2s )8 



(9) 



where E ls = 3 GPa, E 2s = 21.3 MPa, and rj s = 12.4GPa ms. 

In the standard solid model, the initial experimental 
data point is determined by the instantaneous elastic 
modulus E ls . For the indentation that is held for over 
5,000 ms, the indentation force becomes steady at -38 
nN, when the force exerts on the two springs in series. 
In contrast to E ls , E 2s is much smaller, as can be seen 
from the significant force decrease of from -104 to -38 
nN. The tip traveled down 13.2 nm from the beginning 



of indentation. It is noted that for our indentation test, 
the ratio of the maximum indentation depth to the sam- 
ple diameter is less than 10% [48,49]; the substrate effect 
to the elastic modulus calculation is neglected. 

From the determined viscoelastic model, the mechanical 
response of the superlattice under a variety of mechanical 
loads can be predicted. Several simulation results were 
included as follows. 

When the TMV/Ba 2+ superlattice sample undergoes a 
uniformly constant tensile/compressive strain, the stress 
relaxation can be obtained from the standard solid model 
as below 



a(t)=e 0 (E ls +E 2s e- E ^ 



(10) 



where e 0 is the constantly applied strain. 

When the sample undergoes a uniformly constant 
tensile/compressive stress, the strain creep can then be 
obtained as 



e(t) = a 0 



1 

E\s 



1 



-Ei s E 2s t/ rj s (E ls +E 2s ) 



(11) 



where a 0 is the constantly applied stress. 

The stress relaxation vs. applied strains and the strain 
creep vs. applied stresses are shown in Figure 6a,b, re- 
spectively. In Figure 6a, the stress reduces to a steady 
state after ~2 s when the applied strain is -10%. In 
Figure 7b, strain increases to a steady value after ~5 s 
when the applied stress is ~ 1 GPa. 

When the sample is indented with a spherical in- 
denter, the indentation depth history can be analytically 
obtained when a step force is applied. Similar to the 
procedures above where the force history of Equation 
(5) is obtained, a step force function is used as input, 
and the creep indentation depth history function can 
be derived as 

c 3(F 0 + 2ttwR) Gis + G 2s 



d(t) 



A f + —e is + v 

n I* 



2G l5 (l + 6/C l5 ) 

+ G 2s ) : 



Gl s +G 2s )t 



where F 0 is the step force, A f 



4G l4 



■ + 



3/C L 



(12) 



B 



, _ 4GisG 2 2s rj s (l-Gis-G 2 s)-3Ki s rj s GisG 2 s(Gis + 2G 2s + 1) + 3Kisrj s (G 2 2s -G 2 ls -G 



G 2s (Gis + G 2s )(Gis + 2G 2s ) 

_ 4 ! G ls G 2s r} s -3K ls r} s G ls -3K ls r} s Gl s 4G ls G 2s rj s 3K ls rj s G h 



3K ls rj s (G 2s -G ls ) 3K ls rj s G ls 



(Gig + G 25 ) 2 (Gi 5 + 2G 2s ) (Gis + 2G 2s )(Gi s + G 2s ) (G ls + G 2s f G 2s (G ls + G 2s )(G ls + 2G 2s ) G 2s (Gi s + G 2s ) 



2(1 + vis 



3(l-2v ls ) ' 
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The indentation force history has been obtained in 
Equation (5), where the elastic shear modulus G x as a 
combined elastic response of two springs shown in 
Figure 2(b) should be replaced by G ls of one spring only. 
Then, the simulated curves for the two situations can be 
found in Figures 6c,d. It is concluded that the creep 
depth variation under different forces gets larger through 
creep while the indentation force variation under differ- 
ent depths gets smaller through relaxation. Particularly, 
in Figure 6d, the force finally decreases to negative 
values, which represent attractive forces. The attraction 
cannot be found when G l5 and G 2s are very small. This 
phenomenon can be interpreted by the conformability 
of materials determined by the elastic modulus. When 
Gis an d G 2s get smaller, the materials are more con- 
formable. Accordingly, in the final equilibrium state, 
the materials around the indenter tend to be more de- 
formable to enclose the spherical indenter. This will re- 
sult in a smaller attraction. 

In addition, the example of shear dynamic experiment 
is simulated to obtain the storage and loss moduli of 
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TMV/Ba 2+ superlattice. The storage and loss shear mod- 
uli are calculated by [42] 



G f (co) = 0) J G s (t) sin otdt 
o 

G"{a) = 0) J G s (t) cos (otdt 



(13) 



(14) 



where G' and G" are storage and loss moduli, respect- 
ively, o) is the angular velocity which is related to 
the frequency of the dynamic system, and G s (t) = 
Gis + G 2s e~ G2st ^ is the shear stress relaxation modu- 
lus, determined by the ratio of shear stress and con- 
stant shear strain. 

Based on the relation between the transient and dy- 
namic viscoelastic parameters in Equations (13) and 
(14), the storage and loss shear moduli are finally deter- 
mined to be 



G'(a>) 



G» 



Q) 2 G 2s rj 2 s 



(15) 
(16) 



where G 2s = E 2s I 2(1 + v 2s ). 

Figure 7 shows the curves of storage and loss shear 
moduli vs. the angular velocity. The storage shear 
modulus, G ', increases with the increase of angular vel- 
ocity, while the increasing rate of G' decreases and the 
angular velocity of ~2 rad/s is where the increasing rate 
changes most drastically. However, the loss shear 
modulus, G", first increases and then decreases reach- 
ing the maximum value, -3.9 MPa, at the angular vel- 
ocity of -0.7 rad/s. The storage and loss moduli in 
other cases as uniform tensile, compressive, and inden- 
tation experiments can also be obtained. 

Conclusions 

This paper presented a novel method to characterize the 
viscoelasticity of TMV/Ba 2+ superlattice with the AFM- 
based transient indentation. In comparison with previ- 
ous AFM-based dynamic methods for viscoelasticity 
measurement, the proposed experimental protocol is 
able to extract the viscosity and elasticity of the sample. 
Furthermore, the adhesion effect between the AFM tip 
and the sample was included in the indentation model. 
The elastic moduli and viscosity of TMV superlattice 
were determined to be E ls = 2.14 GPa, E 2s = 21.3 MPa, 
and rj s = 12A GPa-ms. From the characterized viscoelas- 
tic parameters, it can be concluded that the TMV/Ba 2+ 
superlattice was quite rigid at the initial contact and 



then experienced a large deformation under a constant 
pressure. Finally, the simulation of the mechanical be- 
havior of TMV/Ba 2+ superlattice under various loading 
cases, including uniform tension/compression and nano- 
indentation, were conducted to predict the mechanical 
response of sample under different loadings. The storage 
and loss shear moduli were also demonstrated to extend 
the applicability of the proposed method. With the char- 
acterized viscoelastic properties of TMV superlattice, 
we are now able to predict the process of tissue regener- 
ation around the superlattice where the time-dependent 
mechanical properties of scaffold interact with the growth 
of tissue. 

Appendix 

Modeling of adhesive contact of viscoelastic bodies 

The functional equation method was employed to de- 
velop a contact mechanics model for indenting a visco- 
elastic material with adhesion. A modified standard solid 
model was used to extract the viscous and elastic param- 
eters of the sample. 

Several adhesive contact models are available, such as 
Johnson-Kendall-Roberts (JKR) model [50], Derjaguin- 
Muller-Toporov (DMT) model [46], etc. [51-53]. De- 
tailed comparisons can be found in reference [54]. As 
the DMT model results in a simpler differential equa- 
tion, it was used in this study for the simulation to solve 
the indentation on an elastic body with adhesion. 

For the DMT model [46], the relation between the in- 
dentation force F and relative approach S, shown in 
Figure 8, can be expressed as 



F + 2jtwR = 8*E*Vr 



(A.l) 




Figure 8 Schematic of contact between a rigid sphere and a 
flat surface (cross-section view). 
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where R is the nominal radius of the two contact spheres 
of Ri and R 2 , given by R = RiR 2 /{Ri + #2); the adhesive 
energy density w is obtained from the pull-off force F c) 
where F c = 3jtwR/2; and the reduced elastic modulus E 
is obtained from the elastic modulus E s and Poisson's 
ratio v s of the sample by E* = 4E S / [3(l-v^)] with the 
assumption that the elastic modulus of the tip is much 
larger than that of the sample. 

In Equation (A.l), E\ which governs the contact de- 
formation behavior, is decided by the sample's mechan- 
ical properties. In the functional equation method [43], 
£* needs to be replaced by its equivalence in the visco- 
elastic system, so that the contact deformation behavior 
can be governed by the viscoelastic properties. To 
achieve it, the elastic/viscoelastic constitutive equations 
are needed. 

As a premise of the functional equation method, 
quasi-static condition is assumed so that the inertial 
forces of deformation can be neglected [43,44]. The gen- 
eral constitutive equations for a linear viscoelastic/ elastic 
system in Cartesian coordinate configuration can be 
written as 



P m a kk = Q m e kk 



(A.2) 
(A.3) 



where 5^, cr kk , and e kk are the deviatoric stress, strain, 
mean stress, and strain, respectively. The linear opera- 
tors K Q d , P™, and Q m can be expressed in the form of 



i=0 



i=0 



pm = Yj>T^ <r = £«r 

i=0 OL i=0 



(A.4a) 



(A.4b) 



where i (i = 0, 1, 2,...) is determined by the viscoelas- 
tic model to be selected, t is time, and pf, qf, p™, 
and q™ are the components related to the materials 
property constants, such as elastic modulus and Pois- 
son's ratio etc. 

For a pure elastic system, the four linear operators are 
reduced to 



P d =p d o, Q? = <ti, P m =Po> Q" 



<7o 



(A.5) 



which, according to the elastic stress-strain relations, are 
correlated as 



4 

Po 



2G 



P% 



pd 



3K 



Q n 

pm 



(A.6) 



where G and K are the shear modulus and bulk modu- 
lus, respectively. 



Combining Equation (A.6) with 
E E 



2(1 + v) 



3(l-2v) 



(A.7) 



the reduced elastic modulus can be expressed by the 
elastic linear operators as 



(A.8) 



^Iqip^pi+pipiq™) 

4 ! (Q d P m Q d + 2P d Q m Q d ) 
~ 3(2Q d P m P d + P d P d Q m ) 

Hence, Equation (A.l) becomes 



(2(/p m P d + p d p d Q m ) [F(t) + 2ttwR] 
= — (Q d p m Q d + 2P d Q m Q d )S^{t) 



(A.9) 



To evolve the elastic solution into a viscoelastic solu- 
tion, the linear operators in the viscoelastic system need 
to be determined. To this end, the standard solid model, 
shown in Figure 2(a), was used to simulate the viscoelas- 
tic behavior of the sample, since both the instantaneous 
and retarded elastic responses can be reflected in this 
model, which well describes the mechanical response of 
most viscoelastic bodies. 

It is customary to assume that the volumetric response 
under the hydrostatic stress is elastic deformation; thus, 
it is uniquely determined by the spring in series [55]. 
Hence, the four linear operators for the standard solid 
model can be expressed as 



pi 



(A.10) 



1,0" 



3K, 



where pf 

G 2 ~ El 



Gj+G 2 

Ki 



4 



2GlG 2 



4 



2G x r\ 



Gi+G 2 ' 11 G±+G 2 ' Wi 2(l+vi) ' 

, Ei, E 2 , V\, and v 2 are the elas- 



■ 2(l+v 2 ) ' " Vi 3(l-2vi) : 

tic modulus and Poisson's ratio of the two elastic com- 
ponents, respectively, shown in Figure 2. 

Plugging Equation (A.10) into Equation (A.9), the rela- 
tion between F(t) and 8{t) can be found. The functional 
differential equation that extends the elastic solution of 
indentation to viscoelastic system is obtained 



Y^Ai—^F^+lnwR] 



i=0 



(A.ll) 



where A 0 = 2q 0 + 3/^, A x = pi{3I<i + 2q 0 ) + {3pJ<i + 2q 1 ) ) 
A 2 = Pi(3pil<i + 2q 1 ), B 0 = q 0 (l + 6 IQ, B 1 = q 0 (p x + 6K&1) + 
qi(6Ki + 1), and B 2 = q x (p x + 6K&1). 
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